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(N . Abstract 

CN ' The Cuba library provides new implementations of four general-purpose multi- 

dimensional integration algorithms: Vegas, Suave, Divonne, and Cuhre. Suave is 
t^J- . a new algorithm, Divonne is a known algorithm to which important details have 

been added, and Vegas and Cuhre are new implementations of existing algorithms 
with only few improvements over the original versions. All four algorithms can in- 
tegrate vector integrands and have very similar Fortran, C/C++, and Mathematica 
! interfaces. 

a: 

Ph! 1 Introduction 
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Many problems in physics (and elsewhere) involve computing an integral, and often enough 
this has to be done numerically, as the analytical result is known only in a limited number 
of cases. In one dimension, the situation is quite satisfactory: standard packages, such as 
Quadpack [1], reliably integrate a broad class of functions in modest CPU time. The 
same is unfortunately not true for multidimensional integrals. 

This paper presents the Cuba library with new implementations of four algorithms for 
multidimensional numerical integration: Vegas, Suave, Divonne, and Cuhre. They have 
a C/C++, Fortran, and Mathematica interface each and are invoked in a very similar 
way, thus making them easily interchangeable, e.g. for comparison purposes. All four can 
integrate vector integrands. Cuhre is a deterministic algorithm, the others use Monte Carlo 
methods. 

Vegas is the simplest of the four. It uses importance sampling for variance reduction, 
but is only in some cases competitive in terms of the number of samples needed to reach a 
prescribed accuracy. Nevertheless, it has a few improvements over the original algorithm 
[21 El and comes in handy for cross-checking the results of other methods. 

Suave is a new algorithm which combines the advantages of two popular methods: 
importance sampling as done by Vegas and subregion sampling in a manner similar to 
Miser By dividing into subregions, Suave manages to a certain extent to get around 
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Vegas' difficulty to adapt its weight function to structures not aligned with the coordinate 
axes. 

Divonne is a further development of the CERNLIB routine D151 5J. Divonne works by 
stratified sampling, where the partitioning of the integration region is aided by methods 
from numerical optimization. A number of improvements have been added to this algo- 
rithm, the most significant being the possibility to supply knowledge about the integrand. 
Narrow peaks in particular are difficult to find without sampling very many points, espe- 
cially in high dimensions. Often the exact or approximate location of such peaks is known 
from analytic considerations, however, and with such hints the desired accuracy can be 
reached with far fewer points. 

Cuhre* employs a cubature rule for subregion estimation in a globally adaptive subdi- 
vision scheme [6J. It is hence a deterministic, not a Monte Carlo method. In each iteration, 
the subregion with the largest error is halved along the axis where the integrand has the 
largest fourth difference. Cuhre is quite powerful in moderate dimensions, and is usually 
the only viable method to obtain high precision, say relative accuracies much below 10~ 3 . 

The new algorithms were coded from scratch in C, which is a compromise of sorts 
between C++ and Fortran 77, combining ease of linking to Fortran code with the avail- 
ability of reasonable memory management. The declarations have been chosen such that 
the routines can be called from Fortran directly. The Mathematica versions are based on 
the same C code and use the MathLink API to communicate with Mathematica. 

2 Vegas 

Vegas is a Monte Carlo algorithm that uses importance sampling as a variance-reduction 
technique. Vegas iteratively builds up a piecewise constant weight function, represented 
on a rectangular grid. Each iteration consists of a sampling step followed by a refinement 
of the grid. The exact details of the algorithm can be found in [21 E] and shall not be 
reproduced here. 

Changes with respect to the original version are: 

• Sobol quasi-random numbers [7j rather than pseudo-random numbers are used by 
default. Empirically, this seems to accelerate convergence quite a bit, most noticeably 
in the early stages of the integration. 

From theoretical considerations it is of course known (see e.g. jHj) that quasi-random 
sequences yield a convergence rate of O ( \og nd n s /n s ), where is the number of 
dimensions and n s the number of samples, which is much better than the usual 
0(l/y/n^) for ordinary Monte Carlo. But these convergence rates are meaningful 
only for large n s and so it came as a pleasant surprise that the gains are considerable 
already at the beginning of the sampling process. It shows that quasi-Monte Carlo 
methods blend well with variance-reduction techniques such as importance sampling. 

*The D from the original name was dropped since the Cuba library uses double precision throughout. 
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Similarly, it was not clear from the outset whether the statistical standard error 
would furnish a suitable error estimate since quasi-random numbers are decidedly 
non-random in a number of respects. Yet also here empirical evidence suggests that 
the standard error works just as well as for pseudo-random numbers. 

• The present implementation allows the number of samples to be increased in each 
iteration. With this one can mimic the strategy of calling Vegas with a small number 
of samples first to 'get the grid right' and then using an alternate entry point to 
perform the 'full job' on the same grid with a larger number of samples. 

• The option to add simple stratified sampling on top of the importance sampling, as 
proposed in the appendix of [2], has not been implemented in the present version. 
Tests with the Vegas version from which contains this feature, showed that con- 
vergence was accelerated only when the original pseudo-random numbers were used 
and that with quasi-random numbers convergence was in fact even slower in some 
cases. 

Vegas' major weakness is that it uses a separable (product) weight function. As a conse- 
quence, Vegas can offer significant improvements only as far as the integrand's characteristic 
regions are aligned with the coordinate axes. 

3 Suave 

Suave (short for subregion- Adaptive VEgas) uses Vegas-like importance sampling combined 
with a globally adaptive subdivision strategy: Until the requested accuracy is reached, 
the region with the largest error at the time is bisected in the dimension in which the 
fluctuations of the integrand are reduced most. The number of new samples in each half 
is prorated for the fluctuation in that half. 

A similar method, known as recursive stratified sampling, is implemented in Miser 
Miser always samples a fixed number of points, however, which is somewhat undesirable 
since it does not stop once the prescribed accuracy is reached. 

Suave first samples the integration region in a Vegas-like step, i.e. using importance 
sampling with a separable weight function. It then slices the integration region in two, as 
Miser would do. Suave does not immediately recurse on those subregions, however, but 
maintains a list of all subregions and selects the region with the largest absolute error for 
the next cycle of sampling and subdivision. That is, Suave uses global error estimation 
and terminates when the requested relative or absolute accuracy is attained. 

The information on the weight function collected in one Vegas step is not lost. Rather, 
the grid from which the weight function is computed is stretched and re-used on the 
subregions. A region which is the result of m — 1 subdivisions thus has had m Vegas 
iterations performed on it. 

The improvements over Vegas and Miser come at a price, which is the amount of 
memory required to hold all the samples. Memory consumption is not really severe on 
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modern hardware, however. The component that scales worst is the one proportional to 
the number of samples, which is 

8(n d + n c + l)n s bytes, 

where is the number of dimensions of the integral, n c the number of components of the 
integrand, and n s the number of samples. For a million samples on a scalar integrand of 
10 variables, this works out to 96 megabytes - not all that enormous these days. 



3.1 Description of the algorithm 

As Suave is a new algorithm, the following description will be fairly detailed. For greater 
notational clarity, n c -dimensional vectors are denoted with a vector arrow (/ ) and n d - 
dimensional vectors with boldface letters (x) in the following, where rid is the dimension 
of the integral and n c the number of components of the integrand. 

The essential inputs are e re \ and e a b s , the relative and absolute accuracies, n^ ew , the 
number of samples added in each iteration, n™ ax , the maximum number of samples allowed, 
and p, a flatness parameter described below. 

Suave has a main loop which calls a Vegas-like sampling step. The main loop is respon- 
sible for subdividing the subregions and maintaining the totals. The sampling step does 
the actual sampling on the subregions and computes the region results. 

3.1.1 Main loop 

1. Initialize the random- number generator and allocate a data structure for the entire 
integration region. Initialize its Vegas grid with equidistant bins. 

2. Sample the entire integration region with n° ew points. This gives an initial estimate 
of the integral I to t, the variance a^ ot , and Xtot- 

3. Find the component c for which r c = cr ctot / max(e abs , e rel / ctot ) is maximal. 
If none of the r c 's exceeds unity, indicate success and return. 

4. If the number of samples spent so far equals or exceeds n™ ax , indicate failure and 
return. 

5. Find the region r with the largest a\. 

6. Find the dimension d which minimizes F c (rf) + F c (r R ), where rf R are the left and 
right halves of r with respect to d. F c {r'l R ) is the fluctuation of the samples that fall 
into r d L R and is computed as 
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where all samples Xj that fall into the respective half are used in the norm/sum and 
the single-sample fluctuation F c is defined as 



F c (x) = w(x) 



/e(x) - I c (r) 



I 



l/c(x)-/ c (r)| 



cr r (r 



This empirical recipe combines the relative deviation from the region mean, (/ — /)/ /, 
with the x value, \f — I\/a, weighted by the Vegas weight w corresponding to sample 
x. Note that the I c and a c values of the entire region r are used. 

Samples strongly contribute to F the more they lie away from the predicted mean 
and the more they lie out of the predicted error band. Tests have shown that large 
values of p are beneficial for 'flat' integrands, whereas small values are preferred if 
the integrand is 'volatile' and has high peaks, p has thus been dubbed a flatness 
parameter. The effect comes from the fact that with increasing p, F becomes more 
and more dominated by 'outliers,' i.e. points with a large F. 

The power 2/3 in Eq. ((!} is also used in Miser, where it is motivated as the exponent 
that gives the best variance reduction ([9 , p. 315). 

7. Refine the grid associated with r, i.e. incorporate the information gathered on the 
integrand in the most recent sample over r into the weight function. This is done 
precisely as in Vegas (see 2J), with the extension that if the integrand has more 
than one component, the marginal densities are computed not from f 2 but from the 
weighted sum' 

c _l c,tot 

8. Bisect r in dimension d: 

Allocate a new region, r^, and copy to Tl those of r's samples falling into the left 
half. Compute the Vegas grid for r L by appropriately "stretching" r's grid, i.e. by 
interpolating all grid points of r with values less than 1/2. 

Set up tr for the right half analogously. 



9. Sample r L with n L = max^ F ( r ^)+F (r R ) n ^ evi ' n ™ n ) an d r R with n R = max(n 
nL,n™ m ) points, where n™ in = 10. 



new 
s 



'It is fairly obvious that scale-invariant quantities must be used in the sum, otherwise the component 
with the largest absolute scale would dominate. It is less clear whether rjo = (J / c dx) 2 = I^ tot , r\\ = 
(/ |/ c | dx) 2 , or rj2 = J dx (or any other) make the best weights. Empirically, rjg turns out to be both 
slightly superior in convergence and easier to compute than 771 and 772 and has thus been chosen in Suave. 

A possible explanation for this is that in cases where there are large compensations within the integral, 
i.e. when J / c dx <C J \ f c \ dx, it is particularly necessary for the overall accuracy that component c be 
sampled accurately, and thus be given more weight in / 2 , and this is better accomplished by dividing / 2 
by the "small" number rjo than by the "large" number 771 or 772- 
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10. To safeguard against underestimated errors, supplement the variances by the differ- 
ence of the integral values in the following way: 



This acts as a penalty for regions whose integral value changes significantly by the 
subdivision and effectively moves them up in the order of regions to be subdivided 
next. 

11. Update the totals: Subtract r's integral, variance, and x 2 -value from the totals and 
add those of r L and r R . 

12. Discard r, put tl and in the list of regions. 

13. Go to Step 3. 

3.1.2 Sampling step 

The function which does the actual sampling is a modified Vegas iteration. It is invoked 
with two arguments: r, the region to be sampled and n m , the number of new samples. 

1. Sample a set of n m new points using the weight function given by the grid associated 
with r. For a region which is the result of m — 1 subdivisions, the list of samples now 
consists of m sets of samples. 

2. For each set of samples, compute the mean Li and variance af. 

3. Compute the results for the region as 






m m m m 



*c = ^2 Wi > cI lc ~ /c W ifihc = ^2 W i,di,c( I i,c ~ h,c) ~ h ^2 W iAkc ~ h,c) ■ 
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4 Divonne 



Divonne uses stratified sampling for variance reduction, that is, it partitions the integration 
region such that all subregions have an approximately equal value of a quantity called the 
spread s, defined as 



where V(r) is the volume of region r. What sets Divonne apart from Suave is that the 
minimum and maximum of the integrand are sought using methods from numerical opti- 
mization. Particularly in high dimensions, the chance that one of the previously sampled 
points lies in or even close to the true extremum is fairly small. 

On the other hand, the numerical minimization is beset with the usual pitfalls, i.e. 
starting from the lowest of a (relatively small) number of sampled points, Divonne will 
move directly into the local minimum closest to the starting point, which may or may not 
be close to the absolute minimum. 

Divonne is a lot more complex than Suave and Vegas but also significantly faster for 
many integrands. For details on the methods used in Divonne please consult the original 
references New features with respect to the CERNLIB version (Divonne 4) are: 

• Integration is possible in dimensions 2 through 33 (not 9 as before). Going to higher 
dimensions is a matter of extending internal tables only. 

• The possibility has been added to specify the location of possible peaks, if such are 
known from analytical considerations. The idea here is to help the integrator find 
the extrema of the integrand, and narrow peaks in particular are a challenge for the 
algorithm. Even if only the approximate location is known, this feature of hinting 
the integrator can easily cut an order of magnitude out of the number of samples 
needed to reach the required accuracy for complicated integrands. The points can be 
specified either statically, by passing a list of points at the invocation, or dynamically, 
through a subroutine called for each subregion. 

• Often the integrand subroutine cannot sample points lying on or very close to the 
integration border. This can be a problem with Divonne which actively searches for 
the extrema of the integrand and homes in on peaks regardless of whether they lie on 
the border. The user may however specify a border region in which integrand values 
are not obtained directly, but extrapolated from two points inside the 'safe' interior. 

• The present algorithm works in three phases, not two as before. Phase 1 performs the 
partitioning as outlined above. From the preliminary results obtained in this phase, 
Divonne estimates the number of samples necessary to reach the desired accuracy 
in phase 2, the final integration phase. Once the phase-2 sample for a particular 
subregion is in, a x 2 test is used to assess whether the two sample averages are 
consistent with each other within their error bounds. Subregions which fail this test 




(2) 
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move on to phase 3, the refinement phase, where they can be subdivided again or 
sampled a third time with more points, depending on the parameters set by the user. 

• For all three phases the user has a selection of methods to obtain the integral estimate: 
a Korobov jTU] or Sobol |7j quasi-random sample of given size, a Mersenne Twister jTTj 
pseudo-random sample of given size, and the cubature rules of Genz and Malik [12J 
of degree 7, 9, 11, and 13 that are also used in Cuhre. The latter are embedded rules 
and hence provide an intrinsic error estimate (that is, an error estimate not based 
on the spread). When this independent error estimate is available, it supersedes the 
spread-based error when computing the total error. Also, regions whose spread-based 
error exceeds the intrinsic error are selected for refinement, too. 

In spite of these novel options, the cubature rules of the original Divonne algorithm 
were not implemented. 

Due to its complexity, the new Divonne implementation was painstakingly tested 
against the CERNLIB routine to make sure it produces the same results before adding 
the new features. 

5 Cuhre 

Cuhre is a deterministic algorithm which uses one of several cubature rules of polynomial 
degree in a globally adaptive subdivision scheme. The subdivision algorithm is similar to 
Suave's (see Sect. I3.1.TJ) and works as follows: 

While the total estimated error exceeds the requested bounds: 

1) choose the region with the largest estimated error, 

2) bisect this region along the axis with the largest fourth difference, 

3) apply the cubature rule to the two subregions, 

4) merge the subregions into the list of regions and update the totals. 

Details on the algorithm and on the cubature rules employed in Cuhre can be found 
in the original references [§]. The present implementation offers only superficial improve- 
ments, such as an interface consistent with the other Cuba routines and a slightly simpler 
invocation, e.g. one does not have to allocate a workspace. 

In moderate dimensions Cuhre is very competitive, particularly if the integrand is 
well approximated by polynomials. As the dimension increases, the number of points 
sampled by the cubature rules rises considerably, however, and by the same token the 
usefulness declines. For the lower dimensions, the actual number of points that are spent 
per invocation of the basic integration rule are listed in the following table. 



number of dimensions 


4 


5 


6 


7 


8 


9 


10 


11 


12 


points in degree-7 rule 


65 


103 


161 


255 


417 


711 


1265 


2335 


4433 


points in degree-9 rule 


153 


273 


453 


717 


1105 


1689 


2605 


4117 


6745 
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6 Download and Compilation 



The source code of the Cuba library can be downloaded as a gzipped tar file from the 
Web site http : / /www . f eynarts . de/ cuba. The archive unpacks into a directory Cuba-1 . 1. 
Change into this directory and type "make" to build the library libcuba. a and the Math- 
Link executables Vegas, Suave, Divonne, and Cuhre. If Mathematica and/or the mcc 
MathLink compiler are not available, type "make lib" to build just the library 

The distribution contains two demonstration programs in Fortran 77 and C, as well as 
the test suite used in Sect. |H1 which is written in Mathematica. 

The code is C99 compliant and compiles flawlessly with the GNU C compiler, versions 
2.95 and higher. Other C compilers may have difficulties with inline functions and variable- 
size arrays, which are C99 extensions. In a pinch, edit the makefile and uncomment the 
line 

CFLAGS += -DNDIM=8 -DNC0MP=2 

This fixes the size of all internal arrays at compile time but, of course, at most 8-dimensional 
integrals of at most 2-component integrands can now be integrated. 

Linking Fortran or C/C++ code that uses one of the algorithms is straightforward, 
just add -lcuba (for the CUBA library) and -lm (for the math library) to the compiler 
command line, as in 

f77 -o myexecutable mysource.f -lcuba -lm 
cc -o myexecutable mysource . c -lcuba -lm 

7 User Manual 
7.1 Usage in Fortran 

Although written in C, the declarations have been chosen such that the routines are directly 
accessible from Fortran, i.e. no wrapper code is needed. In fact, Vegas, Suave, Divonne, 
and Cuhre can be called as if they were Fortran subroutines respectively declared as 



& 
& 
& 



subroutine vegas(ndim, ncomp, integrand, 
epsrel, epsabs, flags, mineval, maxeval, 
nstart, nincrease, 
neval, fail, integral, error, prob) 



& 
& 
& 



subroutine suave(ndim, ncomp, integrand, 
epsrel, epsabs, flags, mineval, maxeval, 
nnew, flatness, 

nregions, neval, fail, integral, error, prob) 
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& 
& 
& 



subroutine divonne(ndim, ncomp, integrand, 
epsrel, epsabs, flags, mineval, maxeval, 
keyl, key2, key3, maxpass, 
border, maxchisq, mindeviation, 



& ngiven, ldxgiven, xgiven, nextra, peakfinder, 
& nregions, neval, fail, integral, error, prob) 

subroutine cuhre(ndim, ncomp, integrand, 
& epsrel, epsabs, flags, mineval, maxeval, 
& key, 

& nregions, neval, fail, integral, error, prob) 

7.1.1 Common Arguments 

• integer ndim (in), the number of dimensions of the integral. 

• integer ncomp (in), the number of components of the integrand. 

• external integrand (in), the integrand. The external subroutine which computes 
the integrand is expected to be declared as 

subroutine integrand (ndim, x, ncomp, f) 

integer ndim, ncomp 

double precision x(ndim) , f (ncomp) 

• double precision epsrel, epsabs (in), the requested relative and absolute accu- 
racies. The integrator tries to find an estimate I for the integral / which for every 
component c fulfills \I C — I c \ ^ max(e abs , e rc \I c ). 

• integer flags (in), flags governing the integration: 

— Bits and 1 encode the verbosity level, i.e. to 3. 

Level does not print any output, level 1 prints 'reasonable' information on the 
progress of the integration, level 2 also echoes the input parameters, and level 
3 further prints the subregion results (if applicable). 

— Bit 2 = 0, all sets of samples collected on a subregion during the various itera- 
tions or phases contribute to the final result. 

Bit 2 = 1, only the last (largest) set of samples is used in the final result. 

— Bit 3 = 0, Sobol quasi-random numbers are used for sampling, 

Bit 3 = 1, Mersenne Twister pseudo-random numbers are used for sampling. 

To select e.g. Sobol quasi-random numbers, last samples only, and verbosity level 2, 
pass 6 = + 4 + 2 for the flags. The higher bits are presently ignored, but should 
be zero for future compatibility. 
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• integer mineval (in), the minimum number of integrand evaluations required. 

• integer maxeval (in), the (approximate) maximum number of integrand evalua- 
tions allowed. 

• integer nregions (out), the actual number of subregions needed (not present in 
Vegas) . 

• integer neval (out), the actual number of integrand evaluations needed. 

• integer fail (out), an error flag: 

— f ail = 0, the desired accuracy was reached, 

— fail = —1, dimension out of range, 

— fail > 0, the accuracy goal was not met within the allowed maximum number of 
integrand evaluations. While Vegas, Suave, and Cuhre simply return 1, Divonne 
can estimate the number of points by which maxeval needs to be increased to 
reach the desired accuracy and returns this value. 

• double precision integral (ncomp) (out), the integral of integrand over the unit 
hypercube. 

• double precision error(ncomp) (out), the presumed absolute error of integral. 

• double precision prob(ncomp) (out), the x 2 -P r obability (not the x 2 -value itself!) 
that error is not a reliable estimate of the true integration error*. 

7.1.2 Vegas-specific Arguments 

• integer nstart (in), the number of integrand evaluations per iteration to start 
with. 

• integer nincrease (in), the increase in the number of integrand evaluations per 
iteration. 

Vegas furthermore allows to store internal parameters for use in subsequent invocations. 
There are two possibilities: 

• It may accelerate convergence to keep the grid accumulated during one integration for 
the next one, if the integrands are reasonably similar to each other. Vegas maintains 
an internal table with space for ten grids for this purpose. The slot in this grid is 
specified by the variable 

*Tb judge the reliability of the result expressed through prob, remember that it is the null hypothesis 
that is tested by the \ 2 test, which is that error is a reliable estimate. In statistics, the null hypothesis 
may be rejected only if prob is fairly close to unity, say prob > .95. 
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integer gridno 

common /vegasgridno/ gridno 

If a grid number between 1 and 10 is selected, the grid is not discarded at the end of 
the integration, but stored in the respective slot of the table for a future invocation. 
The grid is only re-used if the dimension of the subsequent integration is the same 
as the one it originates from. 

• Vegas can also store its entire internal state (i.e. all the information to resume an 
interrupted integration) in an external file. To this end, a file name has to be specified 
in the variable 

character*128 statefile 
common /vegasstate/ statefile 

The state file is updated after every iteration. If, on a subsequent invocation, Vegas 
finds a file of the specified name, it loads the internal state and continues from the 
point it left off. Needless to say, using an existing state file with a different integrand 
generally leads to wrong results. Once the integration finishes successfully, i.e. the 
prescribed accuracy is attained, the state file is removed. 

This feature is useful mainly to define 'check-points' in long-running integrations 
from which the calculation can be restarted. 

7.1.3 Suave-specific Arguments 

• integer nnew (in), the number of new integrand evaluations in each subdivision. 

• double precision flatness (in), the parameter p in Eq. (JTJ), i.e. the type of norm 
used to compute the fluctuation of a sample. This determines how prominently 'out- 
liers,' i.e. individual samples with a large fluctuation, figure in the total fluctuation, 
which in turn determines how a region is split up. As suggested by its name, flatness 
should be chosen large for 'flat' integrands and small for 'volatile' integrands with 
high peaks. Note that since flatness appears in the exponent, one should not use 
too large values (say, no more than a few hundred) lest terms be truncated internally 
to prevent overflow. 

7.1.4 Divonne-specific Arguments 

• integer keyl (in), determines sampling in the partitioning phase: 

keyl = 7, 9, 11, 13 selects the cubature rule of degree keyl. Note that the degree-11 
rule is available only in 3 dimensions, the degree-13 rule only in 2 dimensions. 

For other values of keyl, a quasi-random sample of n\ = |keyl| points is used, where 
the sign of keyl determines the type of sample, 
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— keyl > 0, use a Korobov quasi-random sample, 

— keyl < 0, use a "standard" sample (a Mersenne Twister pseudo-random sample 
if bit 3 of the flags is set, otherwise a Sobol quasi-random sample). 

integer key2 (in), determines sampling in the final integration phase: 

key2 = 7, 9, 11, 13 selects the cubature rule of degree key2. Note that the degree-11 
rule is available only in 3 dimensions, the degree-13 rule only in 2 dimensions. 

For other values of key2, a quasi-random sample is used, where the sign of key2 
determines the type of sample, 

— key2 > 0, use a Korobov quasi-random sample, 

— key2 < 0, use a "standard" sample (see description of keyl above), 

and n 2 = |key2| determines the number of points, 

— n 2 ^ 40, sample n 2 points, 

— n 2 < 40, sample n 2 n nccd points, where n nee d is the number of points needed to 
reach the prescribed accuracy, as estimated by Divonne from the results of the 
partitioning phase. 

integer key3 (in), sets the strategy for the refinement phase: 
key3 = 0, do not treat the subregion any further. 
key3 = 1, split the subregion up once more. 

Otherwise, the subregion is sampled a third time with key3 specifying the sampling 
parameters exactly as key2 above. 

integer maxpass (in), controls the thoroughness of the partitioning phase: The 
partitioning phase terminates when the estimated total number of integrand evalu- 
ations (partitioning plus final integration) does not decrease for maxpass successive 
iterations. 

A decrease in points generally indicates that Divonne discovered new structures of 
the integrand and was able to find a more effective partitioning, maxpass can be 
understood as the number of 'safety' iterations that are performed before the par- 
tition is accepted as final and counting consequently restarts at zero whenever new 
structures are found. 

double precision border (in), the width of the border of the integration region. 
Points falling into this border region will not be sampled directly, but will be extrap- 
olated from two samples from the interior. Use a nonzero border if the integrand 
subroutine cannot produce values directly on the integration boundary. 
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• double precision maxchisq (in), the maximum x 2 value a single subregion is al- 
lowed to have in the final integration phase. Regions which fail this x 2 test and whose 
sample averages differ by more than mindeviation move on to the refinement phase. 

• double precision mindeviation (in), a bound, given as the fraction of the re- 
quested error of the entire integral, which determines whether it is worthwhile fur- 
ther examining a region that failed the x 2 test. Only if the two sampling averages 
obtained for the region differ by more than this bound is the region further treated. 

• integer ngiven (in), the number of points in the xgiven array. 

• integer ldxgiven (in), the leading dimension of xgiven, i.e. the offset between one 
point and the next in memory. 

• double precision xgiven (ldxgiven, ngiven) (in), a list of points where the inte- 
grand might have peaks. Divonne will consider these points when partitioning the 
integration region. The idea here is to help the integrator find the extrema of the in- 
tegrand in the presence of very narrow peaks. Even if only the approximate location 
of such peaks is known, this can considerably speed up convergence. 

• integer nextra (in), the maximum number of extra points the peak-finder subrou- 
tine will return. If nextra is zero, peakf inder is not called and an arbitrary object 
may be passed in its place, e.g. just 0. 

• external peakf inder (in), the peak-finder subroutine. This subroutine is called 
whenever a region is up for subdivision and is supposed to point out possible peaks 
lying in the region, thus acting as the dynamic counterpart of the static list of points 
supplied in xgiven. It is expected to be declared as 

subroutine peakf inder (ndim, b, n, x) 

integer ndim, n 

double precision b (2, ndim) 

double precision x(ldxgiven,n) 

The bounds of the subregion are passed in the array b, where b(l ,d) is the lower and 
b(2,oD the upper bound in dimension d. On entry, n specifies the maximum number 
of points that may be written to x. On exit, n must contain the actual number of 
points in x. 

In contrast to the other algorithms, Divonne passes the integrand one more argument, i.e. 
the integrand subroutine is really declared as 

subroutine integrand (ndim, x, ncomp, f, phase) 

integer ndim, ncomp, phase 

double precision x(ndim) , f (ncomp) 
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The fifth argument, phase, indicates the integration phase: 

• 0, sampling of the points in xgiven, 

• 1, partitioning phase, 

• 2, final integration phase, 

• 3, refinement phase. 

This information might be useful if the integrand takes long to compute and a sufficiently 
accurate approximation of the integrand is available. The actual value of the integral is only 
of minor importance in the partitioning phase, which is instead much more dependent on 
the peak structure of the integrand to find an appropriate tessellation. An approximation 
which reproduces the peak structure while leaving out the fine details might hence be a 
perfectly viable and much faster substitute when phase . It . 2. 

In all other instances, phase can be ignored and it is entirely admissible to declare the 
integrand with only four arguments. 



7.1.5 Cuhre-specific Arguments 

• integer key (in), chooses the basic integration rule: 

key = 7,9, 11, 13 selects the cubature rule of degree key. Note that the degree-11 
rule is available only in 3 dimensions, the degree-13 rule only in 2 dimensions. 

For other values, the default rule is taken, which is the degree-13 rule in 2 dimensions, 
the degree-11 rule in 3 dimensions, and the degree-9 rule otherwise. 



7.2 Usage in C/C++ 

Being written in C, the algorithms can of course be used in C/C++ directly. The decla- 
rations are as follows: 

typedef void (*integrand_t) (const int *, const double [] , 
const int * , double [] ) ; 

void Vegas (const int ndim, const int ncomp, integrand_t integrand, 
const double epsrel, const double epsabs, 
const int flags, const int mineval, const int maxeval, 
const int nstart, const int nincrease, 
int *neval, int *fail, 

double integral [], double error [], double prob[]) 
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void Suave (const int ndim, const int ncomp, integrand_t integrand, 
const double epsrel, const double epsabs, 
const int flags, const int mineval, const int maxeval, 
const int nnew, const double flatness, 
int *nregions, int *neval, int *fail, 
double integral [], double error [] , double prob[]) 

void Divonne (const int ndim, const int ncomp, integrand_t integrand, 
const double epsrel, const double epsabs, 
const int flags, const int mineval, const int maxeval, 
const int keyl, const int key2, const int key3, 
const int maxpass, const double border, 
const double maxchisq, const double mindeviation, 
const int ngiven, const int ldxgiven, double xgiven[] , 
const int nextra, 

void (*peakf inder) (const int *, const double [] , int *, double [] ) , 

int *nregions, int *neval, int *fail, 

double integral [], double error [] , double prob[]) 

void Cuhre (const int ndim, const int ncomp, integrand_t integrand, 
const double epsrel, const double epsabs, 
const int flags, const int mineval, const int maxeval, 
const int key, 

int *nregions, int *neval, int *fail, 

double integral [], double error [] , double prob[]) 

These prototypes are contained in cuba.h which should (in C) or must (in C++) be 
included when using the Cuba routines. The arguments are as in the Fortran case, with the 
obvious translations, e.g. double precision = double. Note, however, the declarations 
of the integrand and peak-finder functions, which expect pointers to integers rather than 
integers. This is required for compatibility with Fortran. 

For convenience, the Divonne prototype glosses over the fact that Divonne passes an 
optional fifth argument to the integrand (see end of Sect. I7.1.4|) . Usually the integrand 
is declared with only four arguments since this extra information is not needed. With 
the 'correct' prototype, the compiler would only generate unnecessary warnings (in C) 
or errors (in C++). In the rare cases where the integrand really has five arguments, an 
explicit typecast to integrand_t must be used in the invocation of Divonne. 

The global variables for the grid number and state file used in Vegas (see Sect. I7.1.2JI 
are also defined in cuba.h as 

extern int vegasgridno_ ; 
extern char vegasstate_ [128] ; 
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7.3 Usage in Mathematica 

The Mathematica versions are based on essentially the same C code and communicate with 
Mathematica via the MathLink API. When building the package, the executables Vegas, 
Suave, Divonne, and Cuhre are compiled for use in Mathematica. In Mathematica one 
first needs to load them with the Install function, as in 

Install ["Divonne"] 

which makes a Mathematica function of the same name available. These functions are used 
almost like NIntegrate, only some options are different. For example, 

Vegas [x~2/ (Cos [x + y + 1] +5), {x,0,5}, {y,0,5}] 

integrates a scalar function, or 

Suave [{Sin [z] Exp[-x~2 - y~2] , 

Cos[z] Exp[-x~2 - y~2]}, {x,-l,l}, {y,-l,3}, {z,0,l}] 

integrates a vector. As is evident, the integration region can be chosen different from the 
unit hypercube. Innermore boundaries may depend on outermore integration variables, 
e.g. Cuhre [1, {x,0,l}, {y,0,x}] gives the area of the unit triangle. 

The functions return a list which contains the results for each component of the inte- 
grand in a sublist {integral estimate, estimated absolute error, % 2 probability}. For the 
Suave example above this would be 

{{1.1216, 0.000991577, 0.0000104605}, 
{2.05246, 0.00146661, 0.00920716}} 

The other parameters are specified via the following options. Default values are given on 
the right-hand sides of the rules. 

7.3.1 Common Options 

• PrecisionGoal -> 3, the number of digits of relative accuracy to seek, that is, 

1 r\— PrecisionGoal 

e re \ — iu 

• AccuracyGoal -> 12, the number of digits of absolute accuracy to seek, that is, 
e a bs = lcr AccuracyGoal . The integrator tries to find an estimate I for the integral I 
which for every component c fulfills | I c — I c \ ^ max(e abs , s re il c ). 

• MinPoints -> 0, the minimum number of integrand evaluations required. 

• MaxPoints -> 50000, the (approximate) maximum number of integrand evaluations 
allowed. 
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• Verbose -> 1, how much information to print on intermediate results, can take 
values from to 3. 

Level does not print any output, level 1 prints 'reasonable' information on the 
progress of the integration, level 2 also echoes the input parameters, and level 3 fur- 
ther prints the subregion results (if applicable). Note that the subregion boundaries 
in the level-3 printout refer to the unit hypercube, i.e. are possibly scaled with respect 
to the integration limits passed to Mathematica. This is because the underlying C 
code, which emits the output, is unaware of any scaling of the integration region, 
which is done entirely in Mathematica. 

• Final -> All, whether only the last (largest) or all sets of samples collected on a 
subregion during the various iterations or phases contribute to the final result. 

• PseudoRandom -> False, whether Mersenne Twister pseudo-random numbers are 
used for sampling instead of Sobol quasi-random numbers. 

• Regions -> False, whether to return the tessellation of the integration region (thus 
not present in Vegas, which does not partition the integration region). 

If Regions -> True is chosen, a two-component list is returned, where the first 
element is the list of regions, and the second element is the integration result as 
described above. Each region is specified in the form Region [x\\ , x UT , res , df] , where 
x\\ and x ur are the multidimensional equivalents of the lower left and upper right 
corner, res is the integration result for the subregion, given in the same form as the 
total result but with the x 2 value instead of the x 2 probability, and df are the degrees 
of freedom corresponding to the x 2 values. 

Cuhre cannot state a x 2 value separately for each region, hence the x 2 values and 
degrees of freedom are omitted from the Region information. 

• Compiled -> True, whether to compile the integrand function before use. Note two 
caveats: 

- The function values still have to pass through the MathLink interface, and in the 
course of this are truncated to machine precision. Not compiling the integrand 
will thus in general not deliver more accurate results. 

— Compilation should be switched off if the compiled integrand shows unexpected 
behaviour. As the Mathematica online help points out, "the number of times 
and the order in which objects are evaluated by Compile may be different from 
ordinary Mathematica code." 

7.3.2 Vegas-specific Options 

• NStart -> 1000, the number of integrand evaluations per iteration to start with. 

• NIncrease -> 500, the increase in the number of integrand evaluations per iteration. 
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• GridNo -> 0, the slot in the internal grid table. 

It may accelerate convergence to keep the grid accumulated during one integration for 
the next one, if the integrands are reasonably similar to each other. Vegas maintains 
an internal table with space for ten grids for this purpose. If a GridNo between 1 
and 10 is chosen, the grid is not discarded at the end of the integration, but stored 
for a future invocation. The grid is only re-used if the dimension of the subsequent 
integration is the same as the one it originates from. 

• StateFile -> " " , the file name for storing the internal state. If a non-empty string is 
given here, Vegas will store its entire internal state (i.e. all the information to resume 
an interrupted integration) in this file after every iteration. If, on a subsequent 
invocation, Vegas finds a file of the specified name, it loads the internal state and 
continues from the point it left off. Needless to say, using an existing state file with 
a different integrand generally leads to wrong results. Once the integration finishes 
successfully, i.e. the prescribed accuracy is attained, the state file is removed. 

This feature is useful mainly to define 'check-points' in long-running integrations 
from which the calculation can be restarted. 

7.3.3 Suave-specific Options 

• Mew -> 1000, the number of new integrand evaluations in each subdivision. 

• Flatness -> 50, the parameter p in Eq. (fTj) . i.e. the type of norm used to compute 
the fluctuation of a sample. This determines how prominently 'outliers,' i.e. individ- 
ual samples with a large fluctuation, figure in the total fluctuation, which in turn 
determines how a region is split up. As suggested by its name, Flatness should be 
chosen large for 'flat' integrands and small for 'volatile' integrands with high peaks. 
Note that since Flatness appears in the exponent, one should not use too large val- 
ues (say, no more than a few hundred) lest terms be truncated internally to prevent 
overflow. 

7.3.4 Divonne-specific Options 

• Keyl -> 47, an integer which governs sampling in the partitioning phase: 

Keyl = 7, 9, 11, 13 selects the cubature rule of degree Keyl. Note that the degree-11 
rule is available only in 3 dimensions, the degree-13 rule only in 2 dimensions. 

For other values of Keyl, a quasi-random sample of ri\ = |Keyl| points is used, where 
the sign of Keyl determines the type of sample, 

— Keyl > 0, use a Korobov quasi-random sample, 

— Keyl < 0, use a "standard" sample (a Mersenne Twister pseudo-random sample 
for PseudoRandom -> True, otherwise a Sobol quasi-random sample). 
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Key2 -> 1, an integer which governs sampling in the final integration phase: 

Key2 = 7, 9, 11, 13 selects the cubature rule of degree Key2. Note that the degree-11 
rule is available only in 3 dimensions, the degree-13 rule only in 2 dimensions. 

For other values of Key2, a quasi-random sample is used, where the sign of Key2 
determines the type of sample, 

— Key2 > 0, use a Korobov quasi-random sample, 

— Key2 < 0, use a "standard" sample (see description of Keyl above), 

and ri2 = |Key2| determines the number of points, 

— n 2 ^ 40, sample n 2 points, 

— n 2 < 40, sample n 2 n nccd points, where n nee d is the number of points needed to 
reach the prescribed accuracy, as estimated by Divonne from the results of the 
partitioning phase. 

Key3 -> 1, an integer which sets the strategy for the refinement phase: 
Key3 = 0, do not treat the subregion any further. 
Key3 = 1, split the subregion up once more. 

Otherwise, the subregion is sampled a third time with Key3 specifying the sampling 
parameters exactly as Key2 above. 

MaxPass -> 5, the number of passes after which the partitioning phase terminates. 
The partitioning phase terminates when the estimated total number of integrand 
evaluations (partitioning plus final integration) does not decrease for MaxPass suc- 
cessive iterations. 

A decrease in points generally indicates that Divonne discovered new structures of 
the integrand and was able to find a more effective partitioning. MaxPass can be 
understood as the number of 'safety' iterations that are performed before the par- 
tition is accepted as final and counting consequently restarts at zero whenever new 
structures are found. 

Border -> 0, the width of the border of the integration region. Points falling into 
this border region are not sampled directly, but are extrapolated from two samples 
from the interior. Use a nonzero Border if the integrand function cannot produce 
values directly on the integration boundary. 

The border width always refers to the unit hypercube, i.e. it is not rescaled if the 
integration region is not the unit hypercube. 

MaxChisq -> 10, the maximum x 2 value a single subregion is allowed to have in the 
final integration phase. Regions which fail this x 2 test and whose sample averages 
differ by more than MinDeviation move on to the refinement phase. 
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• MinDeviation -> .25, a bound, given as the fraction of the requested error of the 
entire integral, which determines whether it is worthwhile further examining a region 
that failed the x 2 test. Only if the two sampling averages obtained for the region 
differ by more than this bound is the region further treated. 

• Given -> O, a list of points where the integrand might have peaks. A point is a list 
of rid real numbers, where rid is the dimension of the integral. 

Divonne will consider these points when partitioning the integration region. The idea 
here is to help the integrator find the extrema of the integrand in the presence of 
very narrow peaks. Even if only the approximate location of such peaks is known, 
this can considerably speed up convergence. 

• NExtra -> 0, the maximum number of points that will be considered in the output 
of the PeakFinder function. 

• PeakFinder -> ({}&), the peak-finder function. This function is called whenever a 
region is up for subdivision and is supposed to point out possible peaks lying in the 
region, thus acting as the dynamic counterpart of the static list of points supplied 
with Given. It is invoked with two arguments, the multidimensional equivalents of 
the lower left and upper right corners of the region being investigated, and must 
return a (possibly empty) list of points. A point is a list of rid real numbers, where 
rid is the dimension of the integral. 

7.3.5 Cuhre-specific Options 

• Key -> 0, chooses the basic integration rule: 

Key = 7,9,11,13 selects the cubature rule of degree Key. Note that the degree-11 
rule is available only in 3 dimensions, the degree-13 rule only in 2 dimensions. 

For other values, the default rule is taken, which is the degree-13 rule in 2 dimensions, 
the degree-11 rule in 3 dimensions, and the degree-9 rule otherwise. 

8 Tests and Comparisons 

Four integration routines may seem three too many, but as the following tests show, all have 
their strengths and weaknesses. Fine-tuning the algorithm parameters can also significantly 
affect performance. 

In the following, the test suite of Genz is used. Rather than testing individual 
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integrands, Genz proposes the following six families of integrands: 



1. Oscillatory: A( x 

2. Product peak: / 2 (x 

3. Corner peak: / 3 (x 

4. Gaussian: / 4 (x 

5. (^-continuous: /s(x 

6. Discontinuous: /g(x 



cos(c • x + 2nwi) 

n d 1 

i=l (Xi - Wi) 1 + c, 



-2 ' 



(1 + C ■ x)™d +1 ' 

exp(— c 2 (x — w) 2 ) , 
exp(— c • |x — w|) , 

for x\ > Wi V x 2 > w 2 

exp(c ■ x) otherwise. 



(3) 



Parameters designated by w are non-affective, they vary e.g. the location of peaks, but 
should in principle not affect the difficulty of the integral. 

Parameters designated by c are affective and in a sense "define" the difficulty of the 
integral, e.g. the width of peaks are of this kind. The q are positive and the difficulty 



increases with ||c||i = Xir=i c «- 

The testing procedure is thus: Choose uniform random numbers from [0, 1) for the Cj 
and Wi. Renormalize c for a given difficulty. Run the algorithms with the integrands thus 
determined. Repeat this procedure 20 times and take the average. 

For comparison, Mathematica's NIntegrate function was included in the test. Unfortu- 
nately, when a maximum number of samples is prescribed, NIntegrate invariably uses non- 
adaptive methods, by default the Halton-Hammersley-Wozniakowski quasi-Monte Carlo 
algorithm. The comparison may thus seem not quite balanced, but this is not entirely true: 
Lacking an upper bound on the number of integrand evaluations, NIntegrate's adaptive 
method in some cases 'locks up' (spends an inordinate amount of time and samples) and 
the user can at most abort a running calculation, but not extract a preliminary result. The 
adaptive method could reasonably be used only for some of the integrand families in the 
test, and it was felt that such a selection should not be done, as the comparisons should 
in the first place give an idea about the average performance of the integration methods, 
without any fine-tuning. 

Table [T] gives the results of the tests as described above. This comparison chart should 
be interpreted with care, however, and serves only as a rough measure of the performance 
of the integration methods. Many integrands appearing in actual calculations bear few 
or no similarities with the integrand families tested here, and neither have the integration 
parameters been tuned to 'get the most' out of each method. 

The Mathematica code of the test suite is included in the downloadable Cuba package. 
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3 


Vegas 


Suave 


Divonne 


Cuhre 


NIntegrate 


1 


162000 ± 


127300 ±32371 


21313 ±11039 


819 ± 


218281 ± 





2 


11750 ± 1795 


13500 ± 1539 


17353 ± 3743 


56238 ±40917 


218281 ± 





3 


16125 ± 2411 


11500 ± 1000 


17208 ± 2517 


1174 ± 444 


218281 ± 





4 


56975 ±11372 


20100 ± 4745 


19636 ± 6159 


22577 ±31424 


218281 ± 





5 


14600 ± 3085 


15250 ± 2337 


21675 ± 4697 


150423 ± 


218281 ± 





6 


19750 ± 4999 


23850 ± 2700 


39694 ±14001 


1884 ± 215 


218281 ± 





n d = 8 


3 


Vegas 


Suave 


Divonne 


Cuhre 


NIntegrate 


1 


153325 ± 20274 


124350 ±35467 


28463 ±31646 


3315 ± 


212939 ±13557 


2 


12650 ± 1987 


21050 ± 4594 


22030 ± 3041 


91826 ±58513 


218281 ± 





3 


24325 ± 3753 


29350 ± 3588 


67104 ±16906 


18785 ±22354 


218281 ± 





4 


38575 ±16169 


29250 ± 8873 


24849 ± 5015 


62322 ±44328 


218281 ± 





5 


15150 ± 2616 


25500 ± 6444 


32885 ± 5945 


151385 ± 


218281 ± 





6 


18875 ± 2512 


40900 ± 7196 


116744 ±32533 


9724 ± 9151 


218281 ± 





n d = 10 


3 


Vegas 


Suave 


Divonne 


Cuhre 


NIntegrate 


1 


156050 ±21549 


129800 ±30595 


32176 ±30424 


7815 ± 


214596 ±164* 


31 


2 


14175 ± 2672 


24800 ± 5464 


25684 ± 7582 


144056 ±25983 


218281 ± 





3 


30275 ± 6296 


51150 ±15608 


139737 ±18505 


109150 ±58224 


218281 ± 





4 


29475 ± 10277 


34050 ± 10200 


27385 ± 8498 


105763 ±49789 


218281 ± 





5 


16150 ± 2791 


31400 ± 7715 


44393 ±18654 


153695 ± 


218281 ± 





6 


22100 ± 3085 


74900 ± 32203 


136508 ±17067 


73200 ±64621 


218281 ± 






Test parameters: 

• number of dimensions: n d = 5, 8, 10, 

• requested relative accuracy: e re \ = 10~ 3 , 

• maximum number of samples: n™ ax = 150000, 



integrand difficulties: 



Integrand family j 


1 


2 


3 


4 


5 


6 


Filli 


6.0 


18.0 


2.2 


15.2 


16.1 


16.4 



Table 1: The number of samples used, averaged from 20 randomly chosen integrands from 
each integrand family j defined in Eq. (jSJ) . Values in the vicinity of n™ ax generally indicate 
failure to converge. NIntegrate seems not to be able to stop at around the limit of 
MaxPoints -> ri™ ax , but always samples considerably more points. 
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9 Summary 



The Cuba library offers a choice of four independent routines for multidimensional numer- 
ical integration: Vegas, Suave, Divonne, and Cuhre. They work by very different methods, 
summarized in the following table (MT = Mersenne Twister): 



Routine Basic integration method 



Algorithm type Variance reduction 



Vegas Sobol quasi-random sample 

or MT pseudo-random sample 

Suave Sobol quasi-random sample 

or MT pseudo-random sample 

Divonne Korobov quasi-random sample 
or Sobol quasi-random sample 
or MT pseudo-random sample 
or cubature rules 

Cuhre cubature rules 



Monte Carlo 
Monte Carlo 

Monte Carlo 
Monte Carlo 

Monte Carlo 
Monte Carlo 
Monte Carlo 
deterministic 



importance sampling 

globally adaptive subdivision 

stratified sampling, 

aided by methods from 
numerical optimization 



deterministic globally adaptive subdivision 



All four have a C/C++, Fortran, and Mathematica interface and can integrate vector 
integrands. Their invocation is very similar, so it is easy to substitute one method by 
another for cross-checking. For further safeguarding, the output is supplemented by a x 2 
probability which quantifies the reliability of the error estimate. 

The source code is available from http://www.feynarts.de/cuba and compiles with 
gcc, the GNU C compiler. The C functions can be called from Fortran directly, so there is 
no need for adapter code. Similarly, linking Fortran code with the library is straightforward 
and requires no extra tools. 

The routines in the Cuba library have all been carefully tested, but it would of course be 
folly to believe they are completely error-free. The author welcomes any kind of feedback, 
in particular bug and performance reports, at hahn@feynarts.de. 
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